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FRACTAL EVALUATION OF TUMOR GROWTH MODELS 


Loretta ICHIM!, Radu DOBRESCU” 


Abstract. The paper presents two applications of fractal analysis techniques in order to 
develop new software instruments for cancer research. Essentially the applications refer 
to two types of tumor growth models, the first in the case of avascular tumors, the second 
in the case of vascularized tumors. The fractal evaluation of the accuracy of tumoral 
growth models is discussed, with care to avoid the conflict between the real complexity of 
the biological process and a reductionist approach for simplest modeling. 
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1. Introduction 


Tumor growth is a most complex process, ultimately dependent on tumor cells 
proliferating and spreading in host tissues. A very important implication of the 
spatial and temporal symmetries of tumors is that certain universal quantities can 
be defined which allow the characterization of the tumor growth dynamics. 


Modeling and simulation of tumor growth in competition with the immune system 
is certainly one of the challenging frontiers of applied mathematical which could 
have a great impact both on the quality of life and development of mathematical 
sciences. The common feature of the above mathematical approach is that the 
equation model living matter and the ability of cells to organize their dynamics 
needs to be an essential feature of these mathematical models. 


Tumor evolution is a most complex process involving many different phenomena. 
Understanding the dynamics of cancer growth is one of the great challenges of 
modern science. Solid tumors develop initially as a single mass of cells. These 
divide more rapidly than the cells around them because of a proliferative 
advantage caused by mutation, and a number of genetic pathways responsible for 
these mutations have been identified over the last decade [1], [2]. 


Because there are three distinct stages (avascular, vascular, and metastatic) to 
cancer development, researchers often concentrate their efforts on answering 
specific questions on each of these stages. Nevertheless, when attempting to 
model any complex system it is wise to try and understand each of the 
components as well as possible before they are all put together. 
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The avascular stage of tumor growth is characterized by small tumors, which gain 
the nutrients and oxygen they need for survival and growth by diffusion from 
external blood vessels. Since there are no blood vessels within the tumor to supply 
the mass needed for such volume expansion, this must also enter through the 
tumor’s periphery. 


An individual tumor cell has the potential, over successive divisions, to develop 
into a cluster of tumor cells. Further growth and proliferation lead to the 
development of an avascular tumor consisting of approximately 10° cells, which 
feed on oxygen and other nutrients present in the local environment. 


Angiogenesis is the process, by which tumors induce blood vessels from the host 
tissue to sprout capillary tips, which migrate towards and ultimately penetrate the 
tumor, providing it with a circulating blood supply and, therefore, an almost 
limitless source of nutrients. 


The vascular growth phase, which follows angiogenesis, is marked by a rapid 
increase in cell proliferation and is usually accompanied by an increase in the 
pressure at the center of the tumor. This may be sufficient to occlude blood 
vessels and, thereby, to restrict drug delivery to the tumor. 


In the earliest stages of development, tumor growth seems to be regulated by 
direct diffusion of nutrients and wastes from and to surrounding tissue. When a 
tumor is very small, every cell receives nourishment by simple diffusion and the 
growth rate is exponential in time. 


However, this stage cannot be sustained because as a nutrient is consumed its 
concentration must decrease towards the center of the tumor. The concentration of 
a vital nutrient at the center will fall below a critical level. 


The main aim of this paper is to develop an accurate growth model for a better 
understanding of process dynamics and the developing of better techniques for the 
prediction of evolution in real instances of cancer. 


2. Related Works 


The process of nutrient consumption and diffusion inside tumors has been 
modeled since the mind 1960°. There have been several reviews ([3], [4]) of this 
area of tumor modeling published over the last few years. However, they all focus 
on different aspects to those we address. 


Most models fall into two categories: A. continuum mathematical models that use 
Space averaging and thus consist of partial differential equations and B. discrete 
cell population models that consider processes on the single cell scale and 
introduce cell-cell interaction using cellular automata type computational 
machinery. 
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A. Continuum cell population models 


Mathematical models describing continuum cell populations and _ their 
development classically consider the interactions between the cell number density 
and one or more chemical species that provide nutrients influence the cell cycle 
events of a tumor cell population. Thus these models typically consist of reaction- 
diffusion equations. One of the best parameterized of these models is due to [5]. 


Early models of nutrient-limited tumor growth calculated the nutrient 
concentration profiles as a function of tumor spheroid radius that was changing 
due to the rate of cell proliferation [6], [7]. The later models have incorporated 
differing degrees of complexity for cell movement. For example, cells can be 
considered to move in either a convective manner or actively in a diffusive 
manner, or in a diffusive/chemotactic manner. 


Most models consider tumor cell proliferation and death to be dependent on only 
one generic nutrient (most often oxygen). The equations describing the 
distribution of molecular species inside the tumor spheroid are classical 
transport/mass conservation equations. 


B. Discrete cell population models 


With the huge advances in biotechnology, large amounts of data on phenomena 
occurring on a single cell scale are now available. This, combined with in vitro 
experiments using tumor spheroids, sandwich culture, etc., and high power 
confocal microscopy that enables tracking of individual cells in space and time, 
has brought about the possibility of modeling single-cell-scale phenomena and 
then using the techniques of up scaling to obtain information about the large-scale 
phenomena of tumor growth. 


There are several up scaling techniques; the most popular ones are cellular 
automata [8], lattice Boltzmann methods [9], agent based [10], extended Potts 
[11] and the stochastic (Markov chain) approach [12]. 


The difficulty with automaton models is realistically modeling cell motion. The 
first step in setting up rules for cell motion is to partition the physical space into 
automaton cells. The simplest partition is to discretize into a regular lattice; 
rectangular lattices are usually chosen for simplicity. 


The second modeling decision is whether the lattice is fixed in time or varies as 
the elements move. It is far simpler to consider a fixed lattice, with each 
automaton cell corresponding to either a biological cell or vacant site, and cells 
able to move into a nearby lattice site containing a vacant site. In particular, while 
the rules of motion for fixed lattices can be formulated simply in terms of cells 
moving between lattice sites, if the lattice is free to move and the cells can grow. 
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A recent three-dimensional cellular automata model, which does not use a regular 
lattice, is that of [13]. The model does not include nutrients or mechanical 
interaction between cells explicitly, but mimics the effect of both in a 
phenomenological way. The authors use random fixed lattice, with the space that 
belongs to a single lattice site consisting of points that are nearer to this site that 
any other lattice site. 


In this model, the proliferation is determined by the distance of the cell from the 
tumor boundary to mimic the effects of nutrient diffusion and consumption; only 
cell within a certain distance from the boundary can proliferate. Similarly, cells a 
certain distance from the boundary become necrotic. 


3. Modeling methods 
A. Application | — modeling avascular tumor growth 
1) General considerations 


Mathematical modeling is an ideal approach for teasing apart mechanisms of 
cancer invasion because it can simultaneously and quantitatively consider 
interactions between multiple factors. In general, simulations may need the use of 
dedicated computer devices, to solve systems, which include biological variables 
in the various transport phenomena related to biological system. The common 
feature of the above mathematical approach is that the equation model living 
matter and the ability of cells to organize their dynamics needs to be an essential 
feature of these mathematical models. 


The biology of tumor micro regions has been investigated experimentally using 
the multicell spheroid model [14], [15]. In experimental setting, the spheroids are 
usually initiated from aggregates consisting of several cells, but as their size 
increases, their growth kinetics become similar to those of tumor in vivo, such as 
micro metastases or pre-vascular primary tumors. The multicell spheroids develop 
layered structure with a central necrotic core surrounded by quiescent cells and a 
tin rim of proliferating cells. Steep gradients in oxygen, glucose and other 
metabolites are also observed in such spheroids. 


The goal of this application was to develop a two-dimensional model to simulate 
the avascular tumor growth based on nutrient consumption. In the same time we 
discuss the fractal analysis as a morphometric measure of the irregular structures 
typical for tumor growth. 


2) The model 


The basic principles included in the model are cell proliferation, quiescent and 
necrosis. Each cell has associated with the velocity, which indicates the direction 
and the distance the cell will move in one time step. 
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There are nine velocity channels in each lattice site: Vo = (0,0), Vi = (,-1), 
V2 = (0,1), V3 = (0,-1), V4 = (1-1), Vs = (-1,0), Ve = (-1,1), V7 = (0,1), Vg = (1,0), 
were Vo is resting channel and Vj, V2, V3, V4, Vs, Vs, V7 and Vg represent moving 
to right, up, left, down and diagonals, respectively. In each lattice site, we allow at 
most one cell (necrotic cells or tumor cells) with each velocity. 


Now in each lattice site one of following reactions can occur at each time step: 


Quiescent: 
: >C,, if and only if (C= 1) (1) 
Ni, Ni; 
Proliferation: 
Ci Cc 
i +N, if and only if (Cij+ Nij <5 and C;,; 21) (2) 
Necrosis: 
fe goog and only if (Ci; > 1) @) 
Ne el : 


where C — tumor cells; N — necrotic cells. 


In order to address the formation of tumor micro regions, we present a two- 
dimensional time-dependent mathematical model in which every tumor cell is 
treated as an individual entity characterized by its own geometry and individually 
controlled cell processes. This model allows one to follow fate of each individual 
cell and to investigate how changes occurring in individual cells can influence 
behavior of the whole tumor tissue. For simplicity, we introduce in our model 
only one external metabolic factor and take explicitly into account only the effect 
of nutrient consumption on cell growth and metabolism. 


Starting from the simplified discrete nutrient equation: 
AN =y7.Cf(N) (4) 


where f(N)=N, if N<1/3 and f(N)=1 otherwise, A, is the second 
numerical difference, C is 1 or 0, depending on whether a cell occupies that lattice 
site or not and 7, is parameter of order 0.1, a scenario of a quasi-random invasion 
process is proposed next. Consider a mnxnn square lattice, with possible values 
nn=16, nn=32, nn=64 or nn=128. The relative nutrient N=n/n,, is the 
invading variable, called invader. Denote (i, j) a place in the lattice, representing 
a cell and consider in (4) the approximation v, =0.3, as specified in [16]. 
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The discrete invader equation (4) can be rewritten, for the cell (i, 7), as: 
N,(K +2) =2N,(k +) —N,(k)+0.3CG, DF(N;(K)) where: (5) 


k is the discrete time, N,, is the invader in the cell @, 7), f(N, (A) =N, (4), if 
N,(k) <1/3, f(N;,(k)) =1 otherwise, CG, j) =1 if the cell (i, j) is occupied by 
the tumour and C(i, 7) =0 otherwise. 


The dynamic model (5) is piecewise linear, with the right hand term depending on 
the values of f and C. For an occupied cell with N,,(k)<1/3, the system (5) is 


linear and free, with the characteristic polynomial: 
UA) =X -2440.7 (6) 
and the eigenvalues 4,, =1+ V3, with A, ols 


For N,,(k) 21/3, the system (5) is forced with a constant term equal to 0.3 in the 
right hand side and has the characteristic polynomial: 


WA HV -241=V +aA+a, (7) 
with the roots 4,, =1. 


Assuming that the nutrient starts from the normal value n,,, i.e. for initial conditions 
in (5) N i (0) = N; (1) =1>1/3, the relative nutrient grows monotonically. 


In this modeling approach, the invasion process has no enemy or partner, so the 
process is not a game. The invasion starts from the center of the lattice, with an 
initial kernel of invaded cells, and at each moment k, if the invader numerical 
gradient reaches, in some already invaded cell, a given limit, then an unoccupied — 
i.e. not belonging to the tumor - or “empty” adjacent cell will be randomly invaded. 


Denote slope a given lower bound for the invader numerical gradient enabling the 
invasion of a neighbor cell. The condition for an event occurrence, i.e. the 
invasion condition from the cell (i, 7) to one of its neighbors (ii, jj) is: 


if{ CG, J) =1} At N, (Kk +2)-N, (kK +1) 2 slope } then . 
{generate rand (i, 7) € (0,1) }A ({jump to (i, jj) = Rule(i, j,rand(i, j)) }, 8) 


mechanism for quasi-random number generation. 


For the cell (i,j) in the nnxnn lattice, with C(i, 7) =1, the pair (ii, jj) can be 
chosen according to the mentioned Rule, so that: 
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(1, J) {GI+D,G-LI.GI-D,G+L DIA 


{C(ii, jj) =0} A{d <ii<nn) ACS jj <nn)}. 0) 


3) Modeling the invasion control as an inhibitor action 


The invasion process is undesired, on one side, and unstable on the other side so it 
is reasonable to try to limit the unwanted invasion by means of an inhibitor, which 
may play the role either of an immune response, as a feedback control, or of a 
treatment. A natural way to build a model of the controlled invasion process is by 
adding a correction variable in the right-hand side of equation (5). Denote M ,(k) 


the current inhibitor for the relative nutrient in a cell (i, 7) . The controlled relative 
nutrient dynamics is given by: 

Nj (k +2) =2N,(k + 1)—N,(k) +0.3CG, J) f(N, (kK) + M i(k) (10) 
The state variables and the control in the model (10) are defined as: 

X(k)=N,(k), x (KVEN (A+), uk) =M,(k). (11) 
The corresponding state realization of the controlled model (10) is: 

X,(k +1) = x,(k) 


X,(k +1) =—x,(k) + 2x,(k) +0.3C7, j) f(x, (k)) + u(k) ) 


Denote x(k) =[x,(k) x,(k)]’ the state of the system (12) at the moment k and 
observe that, similarly to the free model (5), the controlled models (10) and (12) 
are piecewise-linear. In an occupied cell (i,j), 1.e. for CU, 7) =1, if f(x, (k)) =1, 
the dynamics (12) with output y(k) = N,,(k) can be written in the compact form: 


x(k +1) = Ax(k) + b(0.3+ u(k)), y(k) =e" x(k) + du(k) (13) 
with 
a-| a b=]. e"=1 0], d=[0]. (14) 
= es 1 


The characteristic polynomial of the matrix A is yv,(/) = det [al — A] specified in (7). 
The linear system (A,b,c’,d) with the input u(k)=0.3+u(k), k =0,1, 2, ..., is 
controllable so, according to the classical control theory, the closed loop evolution 
can become asymptotically stable by a feedback control u(k)= f x(k) , with 
f’ = Ly, F —€R™ chosen so that the characteristic polynomial 


det [AI-—(A + bf! )] equals some desired polynomial: 
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MW (A) =F + BAF By=(A-AyA-Ano)> Por BeR, (15) 
with io| <1,i=1,2, 
The proposed control problem for the system (12) is slightly different: the objective 


is to bring, in each cell, the relative nutrient to 1=lim, ,,. N(k), for which the 


nutrient equals the “normal” value n =n,, and also to admit that there might appear 
a delay in the inhibitor action. 


The speed of the transition to the desired nutrient value is imposed by the stable 
roots of a given polynomial y, (15). 


B. Application 2 — modeling vascular tumor networks 
1) General considerations 


In the first stages of vascularization, the new vessels are not fully formed. In 
addition, the rate of vascularization may vary significantly between tumors. 


BEE 


K 
; 
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Noes 


{fr 4 
ZZ 


IZ 
Lie 
ZZ. “N\A ASX : 
500um 500um 
arteriovenous network: normal capillary network: carcinoma network: 
Fd= 1.70+0.03 Fd= 1.98+ 0.02 Fd= 1.88+0.04 


Fig. 1. Fractal structure of tumor vessel networks. 
Taking into account that tumor size in the vascular phase of tumor growth is 
directly dependent on the level of vascularization, it was decided to study tumors 
of comparable sizes (approximately 4 mm in diameter). 


Therefore, in tumors the fractal dimension (Fd) was measured for all observed 
vessels taken together. 


In figure 1, we present the real networks and the models based on them — adapted 
from [17]. 


Fractal dimensions were estimated by the box-counting technique [18] where the 
fractal dimension is given by the gradient of the graph of ‘log (box size)’ against 
‘log (number of boxes required to cover the image)’. 
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2) The model 


The computer model of the vascular network starts with a square lattice, invasion 
percolation is implemented by first assigning uniformly distributed random values 
to each point on the lattice. An arbitrary start point is selected as the left bottom 
comer and invades the weakest point the grid at adjacent to the current network 
after each time point; this growth is iterated until the required occupancy is reached. 


The network is then pruned by removing regions with zero flow (dead ends) to 
provide the ‘backbone’, blood vessels are then assumed to connect all adjacent 
lattice sites on this backbone. Selecting the occupancy enables the level of 
geometrical heterogeneity to be controlled as measured by the fractal dimension. 
This model helped to explain why tumor vascular resistance is higher than 
observed in normal vasculature. A foreseeable limitation of this system is fractal 
scaling only applies over a narrow range of length scales and it is not known if the 
scaling extends to larger vessels. 


4. Experimental Results 


For all applications, the models were implemented using Java Development 
Toolkit and the algorithms for computing fractal dimensions were implemented in 
Microsoft Visual C++ 6.0 software. 


A. Application 1 
1) Simulation results 


The simulations conducted on a 100x100 square lattice with central site initially 
defined to contain one cancerous cell. This simulation has been greatly simplified 
by neglecting some effects such as: interaction of healthy cells with cancerous cells, 
the effect of nutrients concentrations and limited volume space for tumor and it 
seems the addition of these effects is not problematic in this simulation (Fig. 2). 


Fig. 2. In vitro growth model versus simulation in two dimensions. 
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Finally, we shall compare the simulated patterns with an in vitro model of tumor 
growth [19] for validation our computational model. In figure 2 we show two 
sections of tumor growth: left panel (a) — in vitro model and right panel (b) — 
simulated pattern. This image is very similar to the patterns exhibited in our 
simulation. At the beginning, we assume that similarities between these in vitro 
growth model and simulated patterns suggest that some of the functional 
properties of cancer cells are similar to those built into our model. 


The size of the lattice is chosen sufficiently large such that the boundaries do not 
influence the tumor growth within the considered time interval. Our model 
estimated fast expansion of tumor cells during the first third of the whole period 
and significant reduction in tumor growth after developing necrotic cell area. 
Finally, the tumor enters into a phase of growth saturation. The percentage of 
proliferating tumor cells is equal to 100% during this time, except of the scattered 
single points that reflect short periods of time, when the newly created daughter 
cells did not yet enter in the new cell cycle. A subpopulation of quiescent cells 
becomes more noticeable at the time when the first necrotic cells arise. After 
subpopulation of necrotic cells arisen the tumor growth is characterized by a fast 
exponential expansion. The percentage of proliferating cells starts to decrease 
with the increasing subpopulations of quiescent and necrotic cells. Three hundred 
and sixty images were automatically generated for each model of tumor growth, 
and their fractal dimensions were estimated using box-counting method. 


2) Computing fractal dimensions 


To measure the fractal dimension and roughness of the tumor boundary we select 
only the cells at the boundary of the tumors. 


1.15 T T T T v 
60 110 160 210 260 310 360 
Time step 


Fig. 3. Plot of the fractal dimensions of the tumor boundary as a function of time steps. 


We define boundary cells as those that have at least one normal neighbor. Fractal 
dimensions were calculated using a box-counting algorithm [18]. See the figure 3 
where the fractal dimensions of the patterns are plotted against the total number of 
time steps. We noticed that the invasive growth is afterwards slowed down while 
the tumor pattern became progressively more compact. This makes us suppose 
that it could be diagnosed in a precocious way the cancer, being based on the 
quantification of the fractal dimension. 
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A. Application 2 
1) Simulation results 


We have developed a computer model based on invasion percolation growth 
process to simulate the geometrical complexity of a tumor vascular network. We 
used the fractal dimension as a quantitative estimator of the spatial complexity of 
the two-dimensional vascular network. Below are some examples of invasion 
percolation clusters in different sizes, with different boundary conditions (Fig. 4). 


Fig. 4. Examples of invasion percolation clusters in different sizes and with different conditions: 
(A) a 100x100 cluster with periodic boundary conditions and Fd = 1.9135; 
(B) a 100x100 cluster with free boundary conditions and Fd = 1.9313; 
(C) a 200x200 cluster with periodic boundary conditions and Fd = 1.9369; 
(D) a 200x200 cluster with free boundary conditions and Fd = 1.9434. 


2) Computing fractal dimensions 


All images were automatically generated and their fractal dimensions were 
estimated. The following table (table 1) includes the average fractal dimension of 
each type of cluster. The average is over 10 runs with different random seeds. All 
of the data are expressed as average fractal dimension + standard deviation. 


Table 1. The average and standard deviation for Fd in different lattice and boundary conditions 


sowiry | ace | arse ec 
Free 100x100 1.92 + 0.02 

Periodic 100x100 1.93 + 0.02 
Free 200x200 1.94 + 0.01 

Periodic 200x200 1.935+0.02 
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5. Conclusions 
We see the role of mathematical modeling in cancer biology as twofold. 


Models can help our intuition, provide a framework for thinking about the 
problem, and make predictions. 


If a model is well parameterized then these predictions can be quantitative 
predictions can be significant. 


1. The proposed abstract invasion model inherits the lattice structure from 
model, but considers only the second-order discrete-time invader equation [20]. 


The simulated evolution of the free invasion process depends strongly on the 
choice of the limit slope. 


On the other side, experiments performed with constant slope value have revealed 
a rather weak evolution variation if the initial conditions are chosen below or 
above the limit 1/3, respectively, for which the characteristic polynomial changes 
from (6) to (7). 


In the proposed controlled model, the correction agent is a growth inhibitor and 
the controlled evolution may reach a stationary state, provided that the inhibition 
is efficient, with the price of an increased complexity. 


The discrete-time approach ensures that for both free and controlled models, the 
numerical simulation is non-blocking and the Zeno phenomenon is avoided. 


In future research, the simple free invasion model can be developed with different 
transition mechanisms, implemented, in (8), by the function Rule. 


Also, simulation of invasion processes in which two or several initial kernels fuse 
gradually, forming a single bigger tumor, may present interest. 


Finally, one can extend the research from this discrete abstraction to a general 
network, thus renouncing to the constraints of a regular lattice. 


2. These results confirm that its whole architecture plays a primary role in 
the quantitative evaluation of the vascular network. 


Therefore, we consider that at the moment to validate our model is necessary 
systematic comparison between the computer simulations and the real data. 


However, both normal and tumor vasculature can more properly be considered 
fractal objects because of their irregular shape (spatial conformation), self-similar 
structure, non-integer dimension and dependence on the scale of observation 
(scaling effect). 


There are many directions in which this work will be taken in the future. 
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In the first research line, we are planning to investigate the transport of diffusible 
substances in tumors, a key process in the growth and treatment of solid tumors. 


For example, adequate oxygen supply is critical for tumor growth but also for 
successful radiation therapy. 


The scale-invariant behavior of vascular networks leads to important insights 
about the transport characteristics of tumors. 


The inherent limitations to oxygen transport in tumor tissues having percolation- 
like vascular networks can be generalized to include most other substances that 
are delivered to target cells by similar convective and diffusive processes. 


Such substances include both nutrients (e.g., glucose) and therapeutic agents (e.g., 
drugs). 


The next line can be, for example, to extend model of invasion percolation in 
three-dimensional lattice. 
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